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ABSTRACT 

The temperature of the intergalactic medium (IGM) is an important factor in deter- 
mining the line- widths of the absorption lines in the Lya forest. We present a method 
to characterise the line-widths distribution using a decomposition of a Lya spectrum 
in terms of discrete wavelets. Such wavelets form an orthogonal basis so the decompo- 
sition is unique. We demonstrate using hydrodynamic simulations that the mean and 
dispersion of the wavelet amplitudes is strongly correlated with both the temperature 
of the absorbing gas and its dependence on the gas density. Since wavelets are also 
localised in space, we are able to analyse the temperature distribution as a function of 
position along the spectrum. We illustrate how this method could be used to identify 
fluctuations in the IGM temperature which might result from late reionization or local 
effects. 

Key words: cosmology: theory - intergalactic medium - hydrodynamics - large-scale 
structure of universe - quasars: absorption lines 



1 INTRODUCTION 

Resonant absorption by neutral hydrogen in the intergalac- 

■ tic medium along the line of sight to a distant quasar is re- 
sponsible for the many absorption lines seen in the Lya for- 

■ est, blueward of the quasar's Lya emission line (Bahcall & 
Salpeter 1965, Gunn & Peterson 1965; see Rauch 1998 for 

, a review). The general properties of these Lya absorption 
lines are remarkably well reproduced by hydrodynamic sim- 
ulations of cold dark matter (CDM) dominated cosmologies 

, (Cen et al. 1994, Zhang, Anninos & Norman 1995, Miralda- 
Escude et al. 1996, Hernquist et al. 1996, Wadsley & Bond 
1996, Zhang et al. 1997, Theuns et al. 1998). 

On large scales where pressure is unimportant, gas 
traces the dark matter and the Lya spectrum can be used 
to infer the underlying density perturbations in the dark 
matter (Croft et al. 1997, Nusser & Haehnelt 1999). On 
small scales however, pressure gradients oppose the infall 
of gas into small potential wells (Jeans smoothing), leaving 
the absorber more extended in space than the underlying 
dark matter. The width of the absorption line is then de- 
termined by residual Hubble expansion across the absorber 
(Hernquist et al. 1996) , Jeans smoothing and thermal broad- 
ening. Theuns, Schaye & Haehnelt (2000) analysed various 
line broadening mechanisms and demonstrated the impor- 
tance of the gas temperature in controlling the line-widths. 

The strong dependence of the small-scale properties of 
the Lya forest on the temperature of the gas allows one 
to reconstruct the thermal evolution of the IGM. The gas 
temperature is set by the balance between adiabatic cool- 
ing caused by expansion and photo-heating by the UV- 



background. This introduces a tight relation between density 
and temperature, T = T (p/ '(p)) 1 ' 1 (Hui & Gnedin 1997). 
The parameters Tb and 7 of this 'equation of state' are very 
sensitive to the reionization history of the IGM (Haehnelt 
& Steinmetz 1998). This is because thermal time scales are 
long in the low density IGM probed by the Lya forest, hence 
that gas retains a memory about the past history of the ion- 
ising background. Consequently, the Lya forest provides us 
with a fossil record of the history of reionization, which can 
be explored by unravelling its thermal history as deduced 
from the Lya forest. 

Schaye et al. (1999, see also Ricotti, Gnedin & Shull 
2000) developed and tested a method to infer To and 7 based 
on the line-widths of the absorption lines. Applying this 
method to high resolution QSO spectra for a range of red- 
shifts, they found (Schaye et al. 2000) that the temperature 
To decreases with decreasing redshift as expected, however, 
there is a large increase in To round z — 3, together with a 
decrease in the value of 7. They attributed this change in the 
equation of state to late reionization of helium II. They also 
noted that the temperature at higher redshifts is still fairly 
high, which might be an indication that we are approaching 
the epoch of hydrogen reionization. 

The method of Schaye et al. to characterise line-widths 
is based on Voigt profile fitting of absorption lines (Webb 
1987, Carswell et al. 1987). The rationale behind fitting ab- 
sorption lines with a Voigt profile is partly historical, and 
stems from earlier theoretical models in which the forest was 
produced by a set of Lya 'clouds' . The line-width of these 
absorbers was assumed to be set by thermal and 'turbu- 
lent' broadening, which would produce a Voigt profile, and 
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line blending was responsible for the lines with large devi- 
ations from the Voigt profile. In the new paradigm of the 
Lya forest absorption in the general IGM is responsible for 
lines, and there is no a priori reason to expect lines to have 
the Voigt shape. 

In this paper we discuss a different method of character- 
ising line- widths, based on discrete wavelets (see e.g. Press et 
al. 1992 for an introduction and further references). Wavelets 
provide an orthogonal basis for a unique decomposition of 
a signal (the spectrum) in terms of localised functions with 
a finite bandwidth. Thus they are a compromise between 
characterising a signal in terms of its individual pixel val- 
ues and in terms of Fourier modes. In the first case, the 
characterisation has no information on correlations between 
different pixels (no frequency information) but perfect po- 
sitional information. A Fourier decomposition, on the other 
hand, has perfect frequency information but no positional 
information. The analysis of a spectrum in terms of wavelets 
has the advantage that one can study the clustering of lines 
('positional information'), as a function of their widths ('fre- 
quency information'). 

The usage of wavelets to analyse QSO spectra was pio- 
neered by Pando & Fang (1996, 1998), who used a wavelet 
analysis of Lya absorption lines to describe the clustering of 
those lines. The wavelet analysis detected large scale struc- 
ture in the Lya forest, which had proved difficult using more 
traditional methods. In contrast to Pando & Fang, we will 
use wavelets to analyse the absorption spectrum directly, 
thereby eliminating the somewhat subjective step of first 
decomposing the continuous spectrum in absorption lines. 
The advantage of this new method is that it allows us to 
objectively characterise the typical width of absorption fea- 
tures as a function of position along the spectrum 

We will show using hydrodynamic simulations that the 
probability distribution of wavelet amplitudes can be used to 
characterise the equation of state of the absorbing medium, 
in terms of the temperature at the mean density, To , and the 
slope, 7, of the temperature-density relation. In addition we 
use the fact that wavelets are localised in position along the 
spectrum, thereby allowing us to detect spatial variations 
in To and/or 7, which might be present as a result of late 
helium II reionization or local effects. 

This paper is organised as follows. In Section ^ we first 
give a brief description of the generation of mock spec- 
tra from our simulations and illustrate the decomposition 
of the spectra in discrete wavelets. The statistics of the 
wavelet amplitudes for different simulations is discussed in 
Section ^ and the results are summarised in Section ^. Re- 
cently, Meiksin (2000) discussed indepently the application 
of wavelets to QSO spectra. 



2 WAVELET ANALYSIS OF MOCK SPECTRA 



2.1 Mock spectra 

We use the LI simulation described before in Theuns et 
al. (2000). Briefly, this is a simulation of a flat, vacuum en- 
ergy dominated cold dark matter model with matter density 
TO = 0.3, baryon fraction Qbh 2 = 0.019 and Hubble con- 
stant Ho = 65 km s _1 Mpc . Density fluctuations in this 
model are normalised to the abundance of galaxy clusters 
(Eke et al. 1996) and we have used CMBFAST (Seljak & 
Zaldarriaga 1996) to compute the appropriate linear transfer 
function. The IGM in this model is photo-ionised and photo- 
heated by the UV-background from QSOs, as computed by 
Haardt & Madau (1996). 

We simulated this cosmological model with a modi- 
fied version of the HYDRA simulation code (Couchman et 
al. 1995), which combines hierarchical P3M gravity (Couch- 
man 1991) with smoothed particle hydrodynamics (SPH, 
Lucy 1977, Gingold & Monaghan 1977). We simulate a pe- 
riodic, cubic box of size 7.7 co-moving Mpc using 128 3 par- 
ticles of each species, which gives us sufficient resolution to 
compute line- widths reliably (Theuns et al. 1998). To in- 
vestigate other effects, we also make use of simulations of a 
model with the same numerical resolution, cosmology and 
thermal history, but with a smaller box size (3.8 Mpc), and 
a set of simulations with a smaller normalisation as = 0.775 
and as = 0.4. 

In the analysis stage, we impose a particular equa- 
tion of state on the gas at low overdensities (p/{p) < 20) 
of the form T = To(p/(p)) 7_1 , varying the values of To 
and 7. We then compute mock spectra that mimick the 
actual observed HIRES spectrum of the z em = 3.0 QSO 
1107+485, discussed by Rauch et al. (1997), using the fol- 
lowing procedure. We divide the observed spectrum in three 
redshifts bins, z = 2.5 - 2.625, 2 = 2.625 - 2.875 and 
z = 2.875 — 3 and scale the mean absorption of the sim- 
ulations at z = 2.5, 2 = 2.75 and z = 3 to the corresponding 
observed value. The simulated spectra are resampled to the 
observed resolution, and convolved with a Gaussian to mim- 
ick instrumental broadening. We have analysed the noise 
statistics of the QSO 1107 spectrum as a function of flux, 
and add noise with these properties to the simulated spec- 
tra. By randomly combining individual sight lines through 
the simulation volume, we generate a single long spectrum 
of length 37 492 kms" 1 . Velocity v is related to redshift z 
via v = c [log e (l + z) — log e (l + Zi)], where c is the speed of 
light, z is redshift and 21 is the redshift where Lya starts to 
be confused with Ly/3 for QSO 1107. In order to perform the 
wavelet analysis, we resample the spectrum to 2 15 =32768 
pixels, equally spaced in velocity. In what follows, we will 
refer to a simulation with a particular equation of state by 
giving To/10 4 K and 7, so the model (1.5,5/3) has the im- 
posed equation of state T = 1.5 x 10 4 (p/(p)) 2/3 . We will 
present results for four equations of state, using To = 1.5 
and 2.2 x 10 4 K and 7 = 1 and 5/3. 



* We will usually refer to absorption features as 'lines', but this 
is just a convenient name for what the eye picks out. The wavelet 
decomposition itself is unique and has no prejudice as to what 
should be considered a line. 



2.2 Wavelets 

The decomposition of a mock spectrum in terms of dis- 
crete wavelets is unique, once a particular wavelet basis has 
been chosen. Here we will use the Daubechies 20 wavelet 
(Daubechies 1988; see e.g. Press et al. 1992 for a general dis- 
cussion on wavelets, and an example of the Daubechies 20 
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Figure 1. Example of a Daubcchics 20 wavelet decomposition of a mock Lya spectrum at z ~ 3. Panel (a): Flux F as a function of 
velocity v for a mock spectrum of QSO 1107. Panel (b): decomposition of F in terms of wavelets with resolutions 2 l ~ 15 X V for i = 9 • • ■ 12 
(from 18.3 to 146.4 km s _1 ) . Panel (c): individual wavelets that make up the curve in (b), for i = 9 (top curve) to i = 12 (bottom 
curve), off-set vertically for clarity. The resolution corresponding to each wavelet is indicated on the right axis. Most lines are detected 
in all shown wavelet resolutions, but only narrow lines are strongly detected at the highest resolution i = 9. 
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Figure 2. Wavelet decomposition of a simulated spectrum at z ~ 3 (panel a) into the wavelet with resolution i = 9 (18.3 km s ), whose 
amplitude A is shown in panel (b). The rms amplitude {A(9, 1000) 2 ), box -car smoothed over 1000 km s is shown in panel (c) (full line). 
The simulated spectrum was made by combining mock spectra from two models with different values of Tq (1.5 X 10 4 and 2.2 X 10 4 K) 
but the same value of 7 = 5/3, in stretches of length 6000 km s — 1 The temperature of this mixed model is shown as the dashed line in 
panel (c) (right axis). There is a strong correlation between the rms wavelet amplitude and the temperature of the absorbing gas, with 
(A 2 ) on average much larger for the cold parts of the spectrum, where To = 1.5 X 10 4 K, than in the hotter parts where To = 2.2 X 10 4 K. 
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wavelet). Just as fast Fourier transforms, (discrete) wavelets 
come in powers of two, but unlike Fourier modes, a given 
wavelet has finite bandwidth and hence corresponds to a 
range of frequencies. Nevertheless we will refer to a wavelet 
of a particular 'resolution', for example quoting its full width 
at half maximum. The simulated spectrum has a length of 
V =37492 km s _1 and the wavelet resolutions correspond 
to 2' -15 x V. Here we will use the exponent i to refer to 
wavelets of a particular resolution, e.g. i = 9 corresponds to 
a wavelet of width 18.3 km s _1 . Analysing a signal in terms 
of the amplitudes of wavelets with different resolutions was 
pioneered in a different context by Mallat (1989). 

An example of a wavelet decomposition of a simulated 
spectrum is shown in Figure [[]. Using wavelets with only 
four resolutions (i = 9 — 12) already gives a relatively good 
description of the strong absorption features in the spec- 
trum. Note how every line in the top panel is 'detected' on 
most resolution levels, indicating that each individual ab- 
sorption line is also made-up of a range of frequencies. This 
is of course because these lines are relatively well approx- 
imated by Voigt profiles, which also have extended band- 
width. However, some lines are only weakly detected in the 
i = 9 narrow wavelet, while some of the narrower lines lead 
to large amplitudes at this high resolution. It is this fea- 
ture, namely that some narrow lines are picked-up strongly 
by the narrow wavelets while the broader lines are not, that 
allows us to characterise objectively the typical line-widths 
of absorption lines. 

For a smaller value To of the IGM temperature, there 
will be a larger fraction of narrow lines in the absorption 
spectrum. For a given pixel at velocity v in the spectrum, 
let 

rv+W/2 

A{v;i,W) = i A{v;i) 2 dv/W (1) 

Jv~W/2 

denote the mean rms amplitude of the wavelet at resolution 
i, box-car smoothed over a window of size W (km s^ 1 ). 
We will usually drop the indices i and v in what follows, 
and assume i = 9 unless stated otherwise. For a spectrum 
with a larger fraction of narrow lines, A will be larger on 
average, hence we can in principle use the statistics of A as 
a measure of To, once the relation between them is calibrated 
with simulations. 

In addition to this mean trend, A will fluctuate along 
the spectrum, due to (random) fluctuations in the strengths 
of lines. Here we give an example showing that averaging A 2 
over a relatively short stretch of spectrum is already enough 
to distinguish between models with different To. This sug- 
gests it might be possible to detect fluctuations in To (and 
7), which might be a relic of a recent epoch of reionization 
or local effects. We will present a more detailed analysis of 
how this can be done below and restrict ourselves here to a 
typical example illustrated in Figure |^. To make the shown 
spectrum, we have combined spectra of the (1.5, 5/3) model 
on scales of 6000 km s~ x with spectra of the 30 per cent 
hotter model (2.2, 5/3), into one long spectrum of length V. 
(In what follows, we will refer to this model as the mixed- 
temperature model.) The rms amplitude A(v; 9, 1000) of the 
i = 9 (18.3 km s _1 ) wavelet, smoothed on 1000 km s , is 
sufficiently different between these two equations of state 
that stretches of the colder model can readily be distin- 
guished from the hotter one as regions with larger A. 
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Figure 3. Cumulative fraction C(< A) of pixels, where the mean 
rms wavelet amplitude A = (A(9, 500) 2 ) of the i = 9 wavelet, box- 
car smoothed over a window of size W = 500 km s — 1 , is less than 
some value, averaged over 100 spectra. The different curves refer 
to different equations of state, as labelled in the figure. Squares 
refer to the mixed-temperature model, obtained from combin- 
ing spectra of model (1.5,5/3) with those of model (2.2,5/3), in 
stretches of length 6000 km s _1 . Models with smaller To and 
shallower equation of state have a larger fraction of pixels with 
large values of A. The mixed-temperature model differs from the 
corresponding single temperature models. 



In this example, both models have been scaled indepen- 
dently to have the same mean optical depth, correspond- 
ing to the observed value for QSO 1107. In reality, regions 
of higher temperature would tend to have smaller optical 
depth because of the T~ ' 7 temperature dependence of the 
recombination coefficient. This would tend to decrease the 
amplitude of the wavelets in the hotter regions even more, 
making it easier to distinguish between hot and cold regions. 



3 WAVELET STATISTICS 

3.1 measuring the equation of state 

In the previous section we showed that the rms amplitude 
of the i = 9 narrow wavelet is strongly anti-correlated with 
the temperature of the absorbing gas. Consequently we can 
characterise the temperature distribution of the IGM over 
the spectrum using the corresponding distribution of wavelet 
amplitudes. For each of 100 realisations of models with a 
specified equation of state, we have computed the cumula- 
tive distribution of A, 

C(< A) = [ P(A)dA, (2) 
Jo 
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Figure 4. Same as Figure 
W = 2000 km s" 1 . 



but for a large smoothing scale 



where P(A) is the probability distribution of A, and we plot 
the mean over 100 realisations, C(< A), in figures ^ and ^ 
for W = 500 and 2000 km s~ , respectively. 

As expected, the colder models are systematically 
shifted to larger values of A, since they contain a large num- 
ber of narrow lines and consequently have larger values of 
A. Note, however, that the dependence on the slope 7 is also 
quite strong, but this may be partly a consequence of using 
the mean density as the pivot point around which we change 
the slope. We have also superposed the mixed-temperature 
model, which stays close to the hot component for small 
values of A before veering away to the locus of the cold 
component for large values of the amplitude. 

Having shown that the mean cumulative distribution 
C(A) depends on the equation of state, we now want to 
characterise how well different models can be distinguished 
from each other, based on a single spectrum. Hence, we want 
to characterise to what extent the cumulative distribution 
Cj (A) for a single spectrum of model j differs from the mean, 
d, for model i. To this end, we compute the dispersion 



(Ci(A)-C 3 (A)) 2 dA. 



(3) 



For a single realisation of a spectrum of model j, Oy is just a 
number. In order to be able to distinguish between two mod- 
els i and j based on a single spectrum, it is necessary that 
the dispersion cr% be much smaller than the mean difference 
crfj between the models. 

Figure |E| shows the cumulative probability distribution 
C(> crfj) for W — 500 for three different equations of state. 
The confidence level at which a single spectrum of the model 
with equation of state say (1.5,5/3) (model j) can be distin- 
guished from the model with equation of state (1.5,1) (model 
i) can be directly read-off from this figure. For example, in 



Figure 5. Cumulative probability distribution of the dispersion 
°ii ~ ((^i ~ ^j) 2 ) f° r a smoothing window W = 500 km s — 1 . 
The first model in the caption refers to the model i for which Cj 
is the mean (over 100 spectra) cumulative distribution of A. The 
second model in the caption refers to model j with cumulative 
distribution Cj for a single spectrum, or?, is a measure of the 
extent that model j resembles model i. For j = i, it characterises 
the dispersion of the cumulative distribution of model i around 
its mean. The shown models (2.2,1) and (1.5,5/3) can both be 
distinguished easily from model (1.5,1) since cr?- tends to be > 
0.004 for most realisations of these models, whereas a realisation 
of model (1.5,1) rarely deviates from its mean to such a large 
extent . 



> 95 per cent of cases afj > 0.004 for i / j, whereas in only 
2 per cent of cases, a model which really has the equation of 
state (1.5,1) will differ from the mean of this model to such 
a large extent. 

A more usual statistic to judge whether a single reali- 
sation of a model is drawn from a given probability distri- 
bution is the Kolmogorov-Smirnov test, based on the maxi- 
mum absolute difference dC = max\Ci{A) — Cj(A)\ between 
two cumulative distributions. The cumulative distribution of 
the KS-statistic is shown in figure |(| where we compare it 
for models (1.5,5/3) and (2.2,1), which resemble each other 
most in figured. For 20 (5) per cent of realisations of model 
(1.5,5/3), dC > 0.1 (dC > 0.12), and at this level of contam- 
ination, 60 (40) per cent of models (2.2,1) have dC larger 
than that. 

Finally, figure (j] illustrates how well the mixed- 
temperature model can be distinguished from either the cold 
or the hot model with 7 = 5/3. This model is most likely 
mistaken with the colder single temperature counterpart. In 
70 (25) per cent of cases, the mixed model has ay > 0.004 
(plj > 0.01). This happens for the cold model in only 10 (5) 
per cent of realisations. 
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Figure 6. Kolmogorov-Smirnov test whether a given realisation 
of model j (second model in the legend) is drawn from the dis- 
tribution of model i (first model in the legend) for i = (1.5,5/3) 
and j = i (full line) and j = (2.2, 1) (dashed line). 



Figure 8. Same as figure [ij but for models with different values of 
the normalisation as, indicated in the legend. Models with smaller 
erg resemble more clustered but hotter models, but the effect is 
relatively weak for realistic values of as- 
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3.2 other effects 

Absorption features are broader in models with a smaller 
amplitude of the dark matter fluctuations (Theuns et al. 
2000), thereby resembling more clustered but hotter models. 
This may lead to a degeneracy between To and as (Bryan 
& Machacek 1999; note that Theuns, Schaye & Haehnelt 
(2000) showed that their Voigt profile analysis does not suf- 
fer from such a degeneracy). For the statistic presented here, 
this degeneracy is not very strong, as shown in figure ^. The 
model with as = 0.775 does not differ much from its more 
clustered counterpart with as ~ 0.9. Only for very low levels 
of clustering, erg = 0.4, is the effect important. All models 
have been scaled to a mean effective optical depth of 0.26 at 
a redshift z = 3. 

Finally we have investigated the influence of the small 
box size in our numerical simulations, and the result is shown 
in figure ^. Lack of long wavelength perturbations decreases 
the observed range in A, as expected, but the effect of this 
purely numerical artifact is relatively weak. 
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Figure 7. Same as Figure ^| but for the mixed-temperature 
model and with W = 2000 km s _1 . The mixed-temperature 
model resembles most strongly its cold constituent, yet can be 
distinguished from it at the > 80 per cent confidence level. 



4 CONCLUSIONS 

Clues to the thermal history of the Universe are hidden in 
the small scale structure of the Lycr forest. There are two 
reasons for this. Firstly, the widths of absorption lines are 
very sensitive to the temperature of gas, and secondly, ther- 
mal time scales are long in the low-density IGM that is re- 
sponsible for the Lycv forest. Since the temperature of the 
photo-ionised IGM is determined by the evolution of the ion- 
ising background, unravelling the thermal history will have 
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Figure 9. Cumulative distribution C(< .4), for a smoothing scale 
of 500 km s _1 (full lines) and 2000 km s — 1 (dashed lines), for a 
simulation box of 3.8 Mpc (thin lines) and 7.7Mpc (thick lines), 
but the same numerical resolution. As before, we have plotted 
C(< ,4) averaged over 100 random realizations of the particular 
model. The model shown is (1.5,5/3) and the spectra are scaled 
to a mean effective optical depth of 0.26 at z = 3. The influence 
of missing long wavelength perturbations on the wavelet statistic 
is small. 



the added benefit of putting strong limits on the sources of 
UV light at high redshifts. 

We have presented a new way of analysing the small 
scale structure of the Lya forest, based on the unique decom- 
position of a spectrum in discrete wavelets. We have shown 
that the rms amplitude (A 2 ) of narrow wavelets (18.3 km 
_1 ) correlates strongly with the temperature of the IGM, 
and also depends on the slope of the equation of state. We 
have quantified to what extent different models can be dis- 
tinguished, using statistics of (A 2 ) 

Our mock spectra have been designed to mimick an 
observed spectrum of QSO 1107+485 as much as possible. 
In particular, we have imposed on our simulated spectra the 



same la :ge scale optical depth fluctuations as are observed in 
QSO llUY, making our mock spectra quite realistic, liven so, 
we can still easily distinguish between models that differ in 
temperature by less than 30 per cent. We have quantified the 
dependence of these statistics on numerical artifacts (miss- 
ing long wavelength perturbations due to the smallness of 
our simulation box) and on the amplitude of the dark matter 
fluctuations (erg). 

Wavelets are also localised in space, making it possi- 
ble to study To and 7 as a function of position along the 
spectrum. We characterised the extent to which we can dis- 
tinguish models with a single value of To from a model with 
temperature fluctuations, as might result from late reioniza- 
tion or local effects. 
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